function [d1, d2, d3] = derivative3D(x1, x2, x3, tensor_values)
% Calculate the derivative of tensor_values with respect to x1, x2, x3 at each
% grid point.  x1, x2, and x3 are grid points, and tensor_values is a
% tensor with dimension (dim(x1), dim(x2), dim(x3))
N1 = numel(x1);
N2 = numel(x2);
N3 = numel(x3);

d1 = zeros(N1, N2, N3);  % derivatives at the first dimension
d2 = zeros(N1, N2, N3);
d3 = zeros(N1, N2, N3);

for( i = 1:N1 )
    d1(i,:,:) = derivative2D( x2, x3,  reshape( tensor_values(i,:,:), N2, N3 )  );
end

for( j = 1:N2 )
    d2(:,j,:) = derivative2D( x1, x3,  reshape(tensor_values(:,j,:), N1, N3)  );
end

for( k = 1:N3 )
    d3(:,:,k) = derivative2D( x1, x2,  reshape(tensor_values(:,:,k), N1, N2)  );
end

